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ABSTRACT 


The problem of radiowave propagation over irregular terrain is solved by using the 
standard parabolic equation method. The ground is characterized by an impedance boundary 
condition and a height profile. A tropospheric boundary condition is used to truncate the 
computational domain. This thesis uses a novel approach of casting all the equations in a 
curvilinear coordinate system. The coordinate system is generated in a simple manner using 
the ground profile data. The equations are solved by the finite difference method using the 
Crank-Nicolson scheme. 

Different numerical values for various important parameters (e.g., step size, location 

_ of tropospheric boundary, the region above the tropospheric boundary, etc.) were used, and 
their effect on the accuracy and computing time are discussed. Validation of the numerical 
results with exact and/or experimental results are presented for different terrain profiles. Both 


perfectly electric conducting (PEC) and lossy impedance surfaces are considered. 
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I. INTRODUCTION 


A. BACKGROUND 

The topic of radiowave propagation over irregular terrain is extremely important in 
ground-to-ground as well as in ground-to-air communications used by the Navy, Air Force, 
and the Army. Similarly, the ability to predict radiowave propagation over irregular terrain 
has a significant impact in determining target detectability in a radar system. The physics of 
propagation is affected by ever-changing atmospheric conditions and by complex terrain 
features on the ground. The link reliability in a communication system or target detectability 
in radar can be significantly affected by the so called 'multipath fading’. The path between a 
transmitter and a receiver is often obstructed by natural or man-made obstacles such as hills, 
buildings, atmospheric layers, trees, rain, fog, etc. In the case of atmospheric multipath fading, 
abnormal propagation of electromagnetic waves resulting from super-refraction or sub- 
refraction can result in severe loss of signals. Reflection multipath fading, which is due to 
interference between the direct and the ground reflected waves depends strongly on the 
terrain geometry and ground constants. It is therefore important to designers and operators 
of communications and radar systems to predict the electromagnetic fields due to radiating 
sources in the troposphere and to assess the effects of environment on radiowave 
propagation. 

Any in-depth understanding of these systems requires a knowledge of the physical 
phenomena that governs low-angle propagation; more specifically, to be able to model 
complex fading phenomena due to refraction, reflection, scattering, and diffraction. Numerous 
analytical methods are available for predicting electromagnetic wave propagation such as 
geometric optics, physical optics, normal mode analysis, and combinations of the above. 
However, they have several limitations for predicting propagation in a complex environment. 
The parabolic equation method has been used to predict radiowave propagation in an 
inhomogeneous atmosphere and over flat terrain, and also for predicting radiowave 


propagation over sloping irregularities. 





In this thesis, we will concentrate solely on the effects of irregular terrain. Since the 
propagation path could extend over several thousands of wavelengths, it is important to have 
a method that is efficient numerically. The parabolic equation method is one such method. 
One advantage of a parabolic partial differential equation (PDE) over an elliptic PDE is that 
in the former case, the field at any location can be computed in terms of the field at a previous 
location. However, this would be accurate only when the waves propagate predominantly in 
the forward direction. In deriving a parabolic PDE, it is assumed that the waves are 
predominantly forward traveling. This is approximately met in a typical radio link where the 
received signal is primarily affected by the nature of the path between the transmitting and 
receiving antennas. Of course there are several situations when this is not true. For example, 
when there is a large obstacle behind a receiving antenna, back-scattering from the obstacle 
will affect the received signal. Nonetheless, the parabolic equation method has been used 
successfully in the past to predict propagation in several scenarios, particularly for 


atmospheric multipath fading. 


B. OBJECTIVE 

In this thesis, we adopt the same parabolic PDE as in the previous approach [Ref. 1] 
to predict radiowave propagation over an irregular terrain. A tropospheric boundary condition 
is used to truncate the computational domain at the top, while an impedance boundary 
condition is used on the uneven terrain to characterize the ground. The key feature of this 
thesis is to use a novel approach of casting all the equations in a curvilinear coordinate 
system. A body fitted coordinate system is generated based on the specification of the ground 
profile. This permits accurate modeling of the boundary conditions which is so vital to the 
success of the model. The parabolic equation method is a full-wave method in that it includes 
all aspects of wave propagation such as forward reflection, refraction, diffraction, and surface 
wave propagation. However, as stated above, it ignores back-scattering. Chapter II presents 
the derivation and formulation of the governing partial differential equation for the standard 
parabolic equation method, the impedance boundary condition, and the tropospheric boundary 


condition. Also presented in Chapter II is the transformation of all the partial differential 











equations to a curvilinear coordinate system. Chapter III details the generation of the 
coordinate system and the numerical procedure for solving the parabolic PDE. The 
performance of the numerical solution is examined in Chapter IV. This includes a study on 
the effects of using different numerical values for various important parameters (e.g., step 
size, location of the tropospheric boundary, the region above the tropospheric boundary, etc.) 
on the accuracy of the solution, and validation of the numerical results with exact and/or 
experimental results for different terrain profiles. Recommendations and conclusions are 
presented in Chapter V. Finally, MATLAB computer codes for the numerical implementation 


are presented in the Appendix. 








Ii. FORMULATION 


In this chapter we present the governing partial differential equation for the fields 
together with the required boundary conditions. We present only the final forms and refer the 
reader to [Ref. 2] for more details. Figure 1 shows a Hertzian electric source placed over an 
irregular, lossy terrain. The terrain 1s characterized by its height profile and an impedance 
boundary condition. The impedance of the ground depends on the ground constants (€,€,, 
itt, and o). We wish to solve the fields at a point on/over the ground in the presence of the 
irregularities. We consider only a 2-dimensional situation where the sources, geometry, and 


all fields are z-invariant. 





A : 
Vv unit normal 


(x0, y0) 


Figure 1. An electric source producing fields over an irregular terrain. 





Section A deals with the impedance boundary condition. In Section B, we present the 
theory on the parabolic PDE. In Section C, we present the tropospheric boundary condition 
required to terminate the computational domain. In any partial differential equation, proper 
imposition of boundary conditions is very critical to the final solution. We desire to solve the 
parabolic PDE by an implicit finite difference scheme. Most of the previous work in this area 
used a cartesian mesh because of its simplicity. However, practical geometries seldom 
conform to cartesian coordinates. Some sort of interpolation is needed near the boundary 
when non-cartesian geometries are encountered, as in the present case. This could result in 
a severe loss of accuracy. To better model the boundary conditions, we solve the equations 
in a curvilinear coordinate system generated by treating the lower irregular boundary as one 


coordinate line. In Section D, we cast all the equations in a curvilinear coordinate system. 


A. IMPEDANCE BOUNDARY CONDITION 

An impedance boundary condition (IBC) relates the tangential components of the 
electric and magnetic fields at the interface of two media. If ¥ is a unit normal vector to the 
boundary, and § is a unit tangential vector as shown in Fig. 1, the boundary condition is [Ref. 
3] 


vx (0x E)--n AY x HF, (1) 
where A, = Z,/n, is the surface impedance normalized to the free-space impedance 
n, = Ju Je, and Z, is the actual surface impedance that is dependent on the media constants 
and the incident angles. The surface impedance is determined from the intrinsic impedance of 
the medium by considering plane wave reflections from the interface. Figure 2 shows the 
interface between two media. The complex propagation constants, y,, y,, and the intrinsic 


impedances 1, and n, can be expressed in terms of the media constants. They are 
2 j : 2 
1 J@p (Gwe) ae “Eos (2) 


, ; 2 
Y> = JOU (oO +jwe€e) = -kKipe , (3) 








6. 


I 


H, 
= Ny» nN, = Nn, — , (4) 
122 
= he plane-wave reflection coefficients for the 
ef. 4] 








(6) 











; Z i = nsec, Horizontal Polarization (7) 
. Z, = 17,0080, Vertical Polarization 
Therefore, 
“ ve cos 2. | ‘ m cos "Wy. 
A, = —-{|1]- —— 3 A, = —-|1- “} . (8) 
re He, €.. He... 
For the special case of normal incidence, w, = 90°, 
ars aye |—_, (9) 
€ -jo 
where o_ = ——. In this thesis, we will use a normalized surface impedance given in Eq. (9) 


@ € 
for all the restlis For a 2-Dimensional case with vertical polarization, the impedance 


boundary condition given by Eq. (1) can be simplified as 


oH 
= 2 jk eA ,H. - 0. IBC Vertical Polarization (10) 
Vv : 


Similarly, we have for the horizontal polarization 


oF 
K =) = 0. IBC Horizontal Polarization (11) 
ov Ae 





For a perfectly conducting material, A,’ = 0 and the impedance boundary condition reduces 


OH, 
to = 0 (Vertical Polarization) and E, = 0 (Horizontal Polarization). 
y 


B. STANDARD PARABOLIC PARTIAL DIFFERENTIAL EQUATION 

For a two-dimensional electric source producing fields in a homogeneous region, all 
quantities are independent of the z-coordinate, and propagation takes place in the xy-plane. 
From Maxwell's equations, we have V x E = -jopH, and V x A - jwek + J. The fields 
could be expressed in terms of the z-component of the magnetic field, H,, in the case of 


vertical polarization (7E, fields), and in terms of E, for the horizontal polarization (7M, 








fields). In a source-free environment, the equation satisfied by the magnetic field is 
V - (VH,) + kH, = 0. Vertical Polarization (12) 
Similarly the equation for the electric field is 
V > (VE) + kE. = 0. Horizontal Polarization (13) 
Equations (12) and (13) may be combined as an equation of the form 
Vy + kw = 0, (14) 


where y = H, for TE Polarization and wy = E, for TM Polarization. Equation (14) is the 
standard Helmholtz equation and 1s elliptic in nature; the field at any one point depends on 
field at every other point in a complicated manner. However, for wave propagation problems, 
an approximate answer can be obtained by the use of a parabolic PDE. In this case, the field 
at a particular range depends on the field at previous range points only. Assuming that the 


wave propagates predominantly in the positive x-direction, we write 


Wix.y) = e**u(x,y), (15) 


where u(x,y) is a Slowly varying function of x. We now impose the restriction that 


ju_| « 2k |u| (16) 
2 
(u, = = , UL = ot ) into Eq. (14) and arrive at 
x ox 
ou Sf ou 17 
ax 2k, ay? i 


Equation (17) is the exact form of the narrow angle parabolic PDE approximation. The 
impedance boundary conditions derived in the previous section can also be expressed in terms 


of the 'u' function : 


u - jk (A, + xu = 0, IBC Vertical Polarization (18) 








A 


& 


l 
ui - a| ra +}: = 0, IBC Horizontal Polarization (19): 


where x, = = =-sinO and wu, = = on the irregular terrain. By defining 
Vv v 


-jk (A; - sin®) Vertical Polarization 


“1° ) jk (—— - sin®) Horizontal Polarization * (20) 
Eqs. (18) and (19) can be combined and written as 
ui+ Cu = 0. (21) 


The parabolic PDE given by Eq. (17) is valid for propagation angles close to the horizontal 
(£15° in practice) [Ref. 5]. 


Cc. TROPOSPHERIC BOUNDARY CONDITION 
Our computational domain consists of the region above the lossy ground. To truncate 
the computational domain on the upper side, we consider a point y = y, high enough and 


impose a boundary condition of the form 
uy + Ou = 6 on y = ¥.. (22) 


To derive this, we start with the parabolic equation u,,— 2jku, = Owith a complex k. 
Consider the problem of determining the field at any point in x > x,,, and y > y,, given the 
initial data on X = Xin, U(XinnY) = f0’), and the boundary data on y = y,, u(x,y,) = g(x). For 
analytical simplicity, we assume k = k, - je, € > 0. The lossless case which we are considering 
in the thesis can be considered as the limit of the lossy case as € —~ 0. Using the Fourier sine 


integral, it is shown in [Ref. 2] that 


ae eo 
ue : lag — e ee ce (23) 


P ae ee ini X-T 
au 








10 





The integrals can be evaluated approximately by replacing the derivatives with 
differences. Figure 3 shows data points on the line y = y, and on the portion of x = x,,, which - 


is above y = y,, 


¥~%N 


- Discrete Data 





of 


0 
U,=U x, x, XX. X, 


Figure 3. _ Discrete derivatives for integral evaluations. 


Consider the evaluation of du/dy|,_,, at x - X,,; = X,.», with initial data on the line x = x,,,. Let 
us assume that this initial data is known on a uniform grid y,, = y, + mAy, m= 0, 1, --. The 


derivative in the interval (y,,,, y,,) can be approximated by the forward difference formula 


re] Ue ey 
ofly) ne yey), (24) 
Oy Ay 


where 4, = U(X4iVn)- hen 


1 { Ym fH1-y, 12 ip,) 5{ Ym ~ Yn) | % 
— S(t )e dt = ma ———— |,| — 
yx TY} k 


Hej) - A [epee 





where F(x) - if * 9 HH" Ie is a complex Fresnel integral [Ref. 6]. 


Now, consider the boundary data on the line y = y,. Assume a non-uniform grid x, 


int> 
Xini T Xp Xing + XQ ~~, Xing + Xp.y, = x. On the interval x - x,,, = (X,), X,,), the derivative can be 


approximated as 


m_ .m-l _ 
g(t) = ———__ (26) 
* ne S| 


1 me 
—— g(t) 1-25 um — ym! Pb PI a 
= e————_— t+ Le ,, 
m=] _ i = 
‘ty =X, x Xu Mma t YXpym Xp Xo Xpy 


Substituting Eq. (25) and (27) into Eq. (23), 


Ou ie *p-% Sik wo *p-% 
oy <, R(X yo i 


7 p-l _ 
: - a eg (28) 
7 m-] x4 x + xy xd (Xo x) 
where Xy = AX, + x) > Xy- xX, = A(x,- X,,) 4 “Ax, ,. Equation (28) can be thought 


of as the discrete version of a continuous boundary condition of the form 





ou r(xyu = s(x), (29) 
dy 


where 


12 





r(x) = i- eS (30) 
Y 7 jx - x 
~* vi 5 —— [fo n(x-x,.) 0,-¥9) — n(x-x,.) Om -»)) 


Sik y a Sik Pl . (31) 
“kr 2 faa caren - x, =: - x m(x - x,,) 


D. TRANSFORMATION TO A CURVILINEAR COORDINATE SYSTEM 
The partial differential equation in Eq. (17) and the boundary conditions in Eqs. (21) 
and (28) will now be transformed to a curvilinear coordinate system. Consider the narrow 


angle parabolic PDE given by Eq. (17) 


fe (32) 


together with the tropospheric boundary condition 


O 
a (x Vo) + P(x,)U(X Vo) = 5(X,Vo)> (33) 
and the impedance boundary condition on the irregular boundary 


u,+cu = 0. (34) 


We will cast all of these equations in a curvilinear coordinate system (€, n) shown in Fig. 4. 
This will permit accurate imposition of boundary conditions. Note that the parabolic nature 
of the equation will not change as a result of the transformation. The irregular terrain 
boundary will map into n = 0 and the upper boundary into 7 = N (integer). The entire 
coordinate system is generated by the specification of the terrain geometry. 

Assuming the transformation x = x(€), y = (&,7), and using the various metrics 


needed [Ref. 7], Eq. (32) becomes 
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a an ; (35) 





Boundary Condition 2 
Physical Domain Computational Domain 


Figure 4. —_ Curvilinear coordinate system. 








Letting 
x Yon) 1 Y 
aa | aera (36) 
n Vy ae ed 
b * | 
3° : (37) 
2k 


up = bu + bu (38) 
The normal derivative on a n = constant line is 


14 





2 2 
x + y x + 


7m 7 4 
Vy Very xe + Ye (x,y, : yx,) 


Using 


x 
———"——— ni Ole | (40) 


2 2 . 2 2 
Xe + Ve et Ve 


the boundary condition at the bottom boundary u, + c,u = 0 gets transformed for x, = 0 to 


cos 8 = 


Vn. 
u,- fm Bove Ou, + cyy,cos Ou = 0. (41) 
4 


Similarly, the boundary condition at y = y, gets transformed to 


(Ege N) + rEg Ny g(E,)#(Ege N) = ¥,S(Eg N). | (42) 
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iil, SOLUTION PROCEDURE 


In the previous chapter we presented the underlying equations and the relevant 
boundary conditions for the computation of fields over an irregular terrain. In Section A of 
this chapter, we present details on the generation of the coordinate system. In Section B, we 
adopt the Crank-Nicolson scheme in the transformed domain to obtain the required finite 


difference equations. 


A. GENERATION OF THE CURVILINEAR COORDINATE SYSTEM 

Consider the computational domain between the terrain profile and a horizontal upper 
boundary. We treat the terrain profile as being comprised of piece-wise linear functions and 
consider a horizontal upper boundary as shown in Fig. 5. We wish to generate a coordinate 
system (€,17) such that the lower and the upper boundaries correspond to constant 1 lines. 
The coordinate lines gradually become horizontal as one proceeds from the lower to the 
upper boundary. Specification of the lower and upper boundaries completely determines our 
coordinate system. Without loss of generality, we will choose the increments between the 
mesh points in the (€,1) plane to be unity. Figure 5(a) shows the mesh points in the physical 
domain and Fig. 5(b) in the computational domain. Figure 6 shows corresponding pairs of 
points on the upper and lower boundaries in the physical and computational domains. The (x, 
y) coordinates of an interior point are generated by using linear interpolation. Letting Ax; = 
(x,,,-x,), Ay,= (4, -y), and AE; = (&;,, - &;), the parametric equations for x and y are [Ref. 
2] 





4 matt: ED 
x( ) ~ x, . A E, j7? (43) 
Ay 
n i n 
y(E.n) ( )| aa "hy ED prods (44) 


for—.< & < &.,,0<n<N. 
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Figure 5. Physical and computational grids. 
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Figure 6. Physical and computational domain segments. 


18 





It is possible to generalize this to any ground profile shape of the form f(x) by writing 


n 7 y n 
N N° N Fie(®)) N° (4) 
with x(€) given as above. However, we choose the piece-wise linear form for its ability to 


model discrete sets of data on the ground. At any interior point, &, < — < &.,,, the various 


metrics are obtained as 


Se ioe x = 0 (46) 
FE,” es 
Ay 1 | Ay 
| i i 
= ] —_— a) 2 > = —x = — ——— E oe . 9 = 0. 


At the junction points € = &, or &,,,, we use the central difference formulas 


Xi 7 %; Ax,, x Ax, 


Ee Bea Ae, ™ 


Ay. + Ay 
4 iv] i (49) 


-E y)-f1- — , 
y(b=E,) ( N) At, > AE, 


Note that the analytical expression yields a discontinuous value for these derivatives. 
We use the analytical expressions only to generate grid points and use the central difference 
formulas to arrive at the derivatives with respect to &. In other words, once the grid points 
are generated on the lines AA’, BB’, ... , etc., we assume that the space is smoothly connected 
through the gnid points. In the numerical implementation using the Crank-Nicolson implicit 
scheme [Ref. 8], the metrics are needed at the midpoint with respect to &, 1.e., at € = &,+ 
A&/2, and the interior point formulas are applicable. For a uniform mesh in the computational 


domain, AE;= 1, n = 9, g = 9, 1,2, ..., N. 
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(50) 


(51) 


— (2) 


(53) 


Figure 7 shows an example of the curvilinear mesh for a Gaussian Shaped ridge on flat 


ground. Note that the vertical increment is constant on any vertical line, but varies from line 


to line. 
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Figure 7. Curvilinear mesh for a Gaussian shaped ridge on flat ground. 
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B. NUMERICAL IMPLEMENTATION USING THE CRANK-NICOLSON | 
SCHEME 


Consider the narrow angle parabolic PDE together with the boundary conditions given 
by Eqs. (38), (41), and (42). We would like to implement the equations using a Crank- 
Nicolson implicit scheme [Ref. 8]. We will write different forms for the internal and boundary 
regions. The grid points for an internal region in the computational domain are shown in Fig. 


8. Given the field data on the line p-1, we would like to predict the field on the line p. 





E=p-1 p-%s Pp 


Figure 8. Internal grid points in computational domain. 


Using the notation of = u(p,g) = u(x,y), the various derivatives (assuming AE = An = 1) 
are 


u (p-'2,9) = us? . ue, (54) 
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plo P Pl oP 
u(p- 2,9) = Aluss* a (ui 4+ ual, 
P Py 4 p-l Pl Pp: 
u(P- 729) = Alu ss 2u)+ uot Uo 2u, + ui). 


Substituting into Eq. (43) and rearranging the terms, we get [Ref. 2] 


1f 4, p | p l b, p 
He + b,) ue = (1 + b)u, + Ms, = > aaa 


Letting 


a - %4(b, + 4b”, 
% 
B a (b,), > 

Y = %4(b, - %b,)?”, 


we have for p = 1, 2, --; g=1, 2, +, N-1 


P P Pp p- p-l p-l 
auo,- (d+ Bu + yu, = -au,, - Ud - Bu: — Yuz). 





(61) 


We will extend the applicability of this equation to include g = 0 and g = N to accommodate 


the derivative boundary conditions on the lower and upper boundaries. For g = 0, we have 


1 1 1 
ou? - (1+ Bu” + yul = -aur - (I - Bue - yul. 
For g = N, we have 


P P P p-| p-l p-1 
Guy, - (1 + B)u,, + Yuy, = -auy, - CI - B)u,, me fete 
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(62) 


(63) 








pl p-Ys p 


Figure 9. Grid points at the lower boundary. 


Figure 9 shows the grid points at the lower boundary. We use central differences for the 


boundary condition in Eq. (41) and combine with Eq. (62) to get 





au? - (1+ Bju? = -auP’ - (1 - Bue, (64) 
where 

a’ = @ + ¥, | (65) 

2sin 0 
B’ = B - 2yy,cos alc, - 20) ; (66) 

3 

2sin 0 

B" = B - 2yy, cos ale + = : (67) 


3 


We have also tried to use one-sided forward difference formulas near the terrain, but the 
results were less accurate than with the central difference scheme. Results for the two 


approaches will be compared in the next chapter. 
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q=N-1 





Figure 10. Grid points at the upper boundary. 


Similarly, using central differences for the tropospheric boundary condition in Eq. (42) and 


combining with Eq. (63) results in 





} 1 r] ] 2 
(1 + Aju, - Yur = (1 - A)uy + yup + days?” (68) 
where 
Jk 
x + 8 a 
Py, nAx 
p-} 
Y=yrta@ 
Ax, =%,- Xp Xoy = X,, + ZAx,, = A(x, + x), 
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 Hie-wo[d (7) AE) 
Y mi TX TX, 





% % 
m m-] me 
+ N°: 
Tw omel ( x ‘ z x + xy x ) Tt Ax, 


We augment Eq. (61) for g = 1, 2, -- , N-1 with Eq. (64) and (68) for g = 0 and q=N 
respectively to define the equations for g = 0, 1, 2, --, N. The system of equations so defined 


can be expressed as a matrix equation of the form 


Pp p-l 
XX Ho XX uo 0 
XXX uP | |X XX uP 
XXX | YX XM | 
: : | (69) 
* 
; i i 
MY, 
Xx X u? XX uP days” 


where X denotes a non-zero entry. The tri-diagonal matrix on the left hand side of Eq. (69) 
can be inverted efficiently to yield a solution on line — = €, in terms of the field values on the 
line — = €,,,. Equation (69) can be used to march forward in range starting from initial data 


specified at § = 0 (x =x,,,). 


2 
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IV. COMPARISON OF RESULTS 


Having presented the basic theory in the previous chapters, we will present numerical 
results in this chapter. The overall accuracy of the solution depends on a number of 
parameters such as the step size, location of the tropospheric boundary (i.e., y, in Fig. 1), the 
region above the tropospheric boundary, and the differencing scheme used at the boundaries. 
In this chapter we will make a systematic study on the dependence of the solution on various 
parameters. Section A.1 compares the results of the forward versus central differencing 
schemes for implementing the lower boundary condition. Sections A.2 and A.3 show the 
effects of horizontal and vertical step sizes, respectively. In Section A.4, we present the effect 
of the placement of the tropospheric boundary condition. In the remaining sections, we 


present comparisons with the results available in the literature for specific obstacles. 


A. PARAMETRIC STUDY OF THE FINITE DIFFERENCE SCHEME 


1. Forward-Difference versus Central Difference Scheme for the Impedance 
Boundary Condition 


In the implementation using the Crank-Nicolson implicit scheme, the derivatives at 
any interior point are approximated using the central difference formulas. This approximation 
is also extended to the lower boundary given in Eq. (69). Alternatively, the derivatives at the 
lower boundary could also be approximated using the forward difference formulas. From Eq. 


(43) and using the forward difference formulas, we have 


Yiu? - uf) + (uP! -uPy - {700 oxin 6 (uy? - up”) 


*e 
+ (cyy,cos | = He = 0. (70) 
Equation (70) can be written as 
aur - (1+ Bu? = -auPr + (1 - Bu, (71) 
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; 2sin 0 
B” = y,,cos | Sy - : (72) 


Numerical results using both the central difference and forward difference schemes 
were used to compute magnetic field due to a vertically polarized line source on the surface 
of a PEC plane. The results of both schemes are compared with the exact solution [Ref. 9] 
and plotted in Fig. 11 (magnitude of magnetic field, Hz, on the surface of the PEC plane), 
Fig. 12 (Hz @ Height = 0.14), and Fig. 13 (Hz @ Height = y, = 20A). It is observed that the 
central difference scheme produces a solution that agrees better with the exact solution. The 
solution using forward differences exhibits spurious oscillations which could not be resolved. 


Henceforth we will adopt the central difference scheme for the subsequent test cases shown. 


2 Variation in Horizontal Step Size (With Fixed Vertical Step Size) 

The Crank-Nicolson scheme has the advantage of being valid (i.e., convergent and 
stable) for all finite values of the horizontal step size for a cartesian mesh. To verify the 
validity of the scheme with respect to horizontal step sizes, numerical results were evaluated 
for three different horizontal step sizes (Ax = 0.14, Ax = 0.24, and Ax = 0.54) with a fixed 
vertical step size of Ay = 0.14. The magnitude of the magnetic field, Hz, due to a vertically 
polarized line source placed at (x0 = 0, y0 = 0) on the surface of a PEC plane are computed 
numerically and plotted in Fig. 14 (Hz on the surface of the PEC plane) and Fig. 15 (Hz @ 
Height = y, = 20A). The plots show that the results are convergent for all the three horizontal 
step sizes, and excellent agreement is observed between the numerical results and the exact 
solution for the surface magnetic field. Minor oscillations are observed for magnitude of the 
magnetic field at the upper boundary. However, the average value of the magnetic field, Hz, 
is still close to the exact solution given by [Ref. 9] 

Sk 


, &-£ : 
H (x,y) = 2 | 7 Soe HR}. (73) 


R 


° 
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where f- kx, X,-kx,, R,= kale - x)? +0 - 9), and H,” is the Hankel 
function of the second kind of order 1. 

3. Variation in Vertical Step Size (With Fixed Horizontal Step Size) 

Next, we investigate the behavior of the Crank-Nicolson scheme with respect to 
vertical step size. In this study, the numerical results were evaluated for three different vertical 
step sizes (Ay = 0.14, Ay = 0.2A, and Ay = 0.5A) with a fixed horizontal step size of Ax = 
0.14. The same profile of a vertically polarized line source placed at (x0 = 0, y0 = 0) on the 
surface of a PEC plane is adopted, and the magnitude of the magnetic field, Hz, are computed 
numerically. The results are plotted in Fig. 16 (Hz on the surface of the PEC plane) and Fig. 
17 (Hz @ Height = y, = 20). Again the plots show that the results at the surface are 
convergent for all the three vertical step sizes tested, and excellent agreement is observed 
between the numerical results and the exact solution for the magnetic field, fz. 

Minor oscillations are still observed at the upper boundary. Poor results are obtained 
with Ay = 0.5A. This is possibly due to the inaccuracies of the upper boundary condition. 
Nonetheless, the numerical computation using step size of Ax = Ay = 0.1A produces results 
which are in good agreement with the exact solution. 

4. Variation in the Tropospheric Boundary Condition Parameter 

In the evaluation of du/dy|,_,. at x =X, ,+%j,, we have considered and assumed that 
initial data is known on a uniform grid y, - y,= mAy, m= 0, 1, --. In practice we stop at some 
height y,,, > Yo. To ascertain how high a point should be considered, we examine the cases 
of extending beyond the tropospheric boundary by 5A, 104, and 204, (1.€., Ys. =Yo + 5A, 
10A, and 20A). The simple test profile of the propagation of the magnetic field, 7z, due to 
a vertically polarized line source placed at (x0 = 0, y0 = 0) on the surface of a PEC plane are 
computed numerically. In this study we choose Ax = Ay = 0.14 and the upper boundary, y, 
= 20A. The results are plotted in Fig. 18 (Hz on the surface of the PEC plane) and Fig. 19 (Hz 
@ Height = y, = 20A). The plots show that the results are almost identical for all the three 
cases. To get an accurate representation of the tropospheric boundary condition, we need 


only to consider data on a vertical line which is within 5A from the upper boundary. As such, 


all subsequent test cases will be computed using a tropospheric boundary with y,,, =y,t5A. 
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Figure 11. Surface magnetic field versus distance. 
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Figure 12. Hz versus distance @ height = 0.1 A. 
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Figure 13. Hz versus distance @ height = y, = 20A. 
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Figure 14. Surface magnetic field versus distance. 
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Figure 15. Hz versus distance @ height = y, = 202. 
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Figure 16. Surface magnetic field versus distance. 
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Figure 17. Hz versus distance @ height = y, = 204A. 
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Figure 18. Surface magnetic field versus distance. 
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Figure 19. Hz versus distance @ height = y, = 204. 
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B. PROPAGATION OVER A PEC PLANE | 
In this section, we examine the effect of the height of the upper boundary, y,, on the 
stability and accuracy of the mane results. The magnitude of the magnetic field, Hz, due 
to a vertically polarized line source placed at (x0 = 0, y0 = 0) on the surface of a PEC plane 
are computed numerically. In this study we choose Ax = Ay=0.1A andy, =y, + SA. 
Results are presented for y,= 5A, 10A, 20A, and 30A. The initial data line is at x, = 5A, and 


the magnetic field is marched over a horizontal distance of 10A. 


1. Magnetic Field Variation with Horizontal Distance 

The variation in the magnitude of the magnetic field, Hz, along the horizontal distance 
on the surface of the PEC plane and at a height of 5A is plotted in Fig. 20 and Fig. 21 
respectively. We noticed that the oscillations are predominant when y, is less than 20A. For 


y, > 20A, the results are in good agreement with the exact solution. 


2. Magnetic Field Variation with Vertical Distance 

Figures 22, 23, and 24 depict the results of the variations of the magnitude of the 
magnetic field, Hz, with vertical height at horizontal distances of 6A, 104, and 15A from the 
transmitter. Again, oscillations are present when the height of the upper boundary is below 
204. More importantly, we observed that for a given horizontal distance, the deviations of 
the numerical results from the exact solution become more pronounced at higher points. This 
phenomenon is expected because the parabolic equation used in this thesis is valid for 
propagation angles close to horizontal (+15°). In fact the numerical result 1s fairly consistent 
with this restriction. For example, at a horizontal distance of 154 from the transmitter (i.e., 
10A from the initial data line), the maximum height for which good results are expected is 
between 0 and 1.7634. From Fig. 24, we see that this 1s true, and deviation between the 
numerical results and the exact solution begins at approximately 1.71. The same trend is also 
observed in Fig. 22 and 23. However, the difference between the numerical results and the 
exact solution increases with distance from the transmitter. This could be due to inaccuracies 
of the tropospheric boundary condition. The error in the tropospheric boundary condition 1s 


accumulative, and we would expect to see higher error at larger distances. 
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Figure 20. Surface magnetic field versus distance. 
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Figure 21. Hz versus distance @ height = 5A. 
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Figure 22. Hz versus height @ distance = 6A. 
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Figure 23. Hz versus height @ distance = 104A. 
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Figure 24. Hz versus height @ distance = 15A. 
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C. PROPAGATION OVER A PEC KNIFE-EDGE 


Consider a perfectly conducting knife edge of height 4, as shown in Fig. 25, set ona 


perfectly conducting plane, with the transmitter and receiver both on the ground at distances 


d, and d,. The attenuation factor, AF, at point B can be found analytically and is given by 


[Ref. 10] 


AF - 4 2 [e jru*p du , 
9) v, 
2d 
u(t) = § u, = u(Esh) = h,|—~ . 
id,d, id,d, 


where 











Tx ye \ Rx 


‘ . 





0, 0, 





d, d, 





Figure 25. Perfectly conducting knife edge between the transmitter at A and the receiver 


at B, both of which are on a perfectly conducting ground. 
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A.’ 


(74) 


(75) 








We consider the case where h, = 2A and d, = 102. For this choice 0, = 11.3°. A vertically 
polarized line source is placed at x0 = 0, y0 = 0 and the magnetic field is determined from an 
initial range of 5A to a distance of 304. The height of the upper boundary is y, = 20A and y,,. 
= Vor DK, 

The variation of the surface magnetic field over the horizontal distance is plotted in 
Fig. 26. Variation of the magnetic field with height at a distance d, = 2A is plotted in Fig. 27. 
Comparison of the attenuation factors at d, = 2A, 104, 154, and 20A (6, = 45°, 11°, 7.6°, 
ee respectively) computed numerically and the exact solution using Eq. (74) is given in 


Table 1. 





Table 1. Attenuation factor for propagation over PEC knife edge. 


The rather large oscillations seen near the surface in Fig. 27 could be due to the non- 


applicability of the parabolic equation for high angles (i.e., 6, > 15°). 


D. PROPAGATION OVER A CIRCULAR BOSS IN A PEC PLANE 

Next, we compare the numerical results for propagation over a semi-circular boss on 
a perfectly conducting plane. This problem typifies radiowave propagation over uneven 
terrain. A vertically polarized line source is placed in front of the semi-circular boss of radius 
b = 0.5A and centered at the origin. The transmitter is at x0 = -0.75A, yO = 0.14 from the 
ongin and the initial data line is at -0.6A from the center of the boss. There is good agreement 
between the exact solution given in [Ref. 9] and the numerical results generated here as 


evidenced in Fig. 28. 
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Figure 26. 
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Normalized magnetic field versus height @ d2 = 24 (PEC knife edge). 


48 





H1z/Hz(f,s)| 










—— Numerical Solution 
© Exact Solution [Ref. 8] 
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Figure 28. Magnitude of the normalized surface fields versus horizontal distance for 
a line source placed over a perfectly conducting plane with a semi-circular boss. 
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E. PROPAGATION OVER A LOSSY IMPEDANCE PLANE 

In the previous sections, we have considered only propagation over PEC surfaces. In 
this section, we will investigate the case of propagation over a lossy impedance plane. The 
ground constant chosen are €, = 10, o, = 180 (corresponding to o = 10 mS/m at 1 MHz). 
The line source is placed at x0 = 0, yO = 0.014. The horizontal step size, Ax = 0.1A, and the 
vertical step size, Ay =0.1A. The initial data line is 5A away from the source, and the solution 
is allowed to propagate for 25 wavelengths. The magnetic field on the initial data line is 


generated from the fields on a flat lossy plane [Ref. 11] 


k : Sk re 
H,- -— | cos OH (ky, + cos OH Er, - 724 |A~-2 "| (76) 
4] tk, (A, + sin®,) 





where 


rey -x)+O-y), r,7 ¥@- x) +ye+y)’, 


gee, weg 


ry Ka 


P Ten 4 
ae er 


r 


cos 0, = : 


1 


The results for the surface magnetic field are plotted in Fig. 29. Comparison is made 
with the data points for the exact solution given in Fig. 22, pp. 52 of [Ref. 9]. There is very 


good agreement between the numerical results and the exact solution. 


F. PROPAGATION OVER A LOSSY GAUSSIAN HILL 
Next, we treat the case of propagation over a Gaussian shaped ridge or hill. The 


height h(x) of the Gaussian hill is defined as 
h(x) = He eo (77) 


where // is the maximum height, c is the location of the maximum, and w is a parameter 
controlling the width of the hill. A vertically polarized line source was placed at x0 = 0, y0 = 
0.014 with H= 1 km, c=5 km, and w = 3 km. The initial data line was at x,,, = 2 km. 
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Figure 29. Surface magnetic field versus distance (lossy plane). 
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Figure 30. Magnitude of the normalized surface fields versus horizontal distance for a 
source placed at the origin in front of a Gaussian hill. 
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The ground constants were e, = 10 and o, = 180 (corresponding to o = 10 mS/m at 1 MHz). 
We plot the results of the normalized magnetic field in Fig. 30 together with the results 
generated from WAGNER [Ref: 12]. Although there is minor differences in the magnitude 
of the normalized magnetic field, it is seen that the trend in the two results are in good 
agreement. Both results predict large peak around 3.9 km from the source, which is attributed 
to focusing by the concave surface of the hill. A secondary peak at 4.5 km is also predicted 


by the numerical method. 


G. PROPAGATION OVER A PEC ISOSCELES HILL 

Experimental results for propagation over a PEC isosceles triangular hill are given in 
[Ref. 10]. In [Ref. 10], an irregular path consisting of an isosceles hill of height 1.234 was 
constructed from aluminum sheets whose surface impedance was taken to be zero. The 


dimensions of the sloping side of the hill are shown in Fig. 31. 


Initial Data Line Hill H = 3.85 cm 
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(10.24A) 


Figure 31. Perfectly conducting isosceles triangular hill of height 1.234 and baselength 
20.164. 
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Figure 32. Magnitude of the normalized fields across a perfectly conducting isosceles 
triangular hill. 
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Using the vertical line sources at a distance of 10.24A from the beginning of the hill, 
and a frequency of 9.6 GHz, the numerical solution is computed and marched from an initial 
data line of 5.124. The normalized magnetic field taken at a height 4, = 0.164 above the 
irregular surface is plotted in Fig. 32. We have computed the attenuation factor for the 


experimental values given in Fig. 4, pp. 39 of [Ref. 10] as follows 


| : (0 | ( | - 
AF (dB) = 201 ——j +. 10] — 
(dB) oo( = og : (78) 


ri ri 


and compared them to the results obtained from our numerical method. In our case, we have 
chosen r; = 5.12A to coincide with the initial data line. It can be seen that there is excellent 


agreement between the experimental data and the results from the numerical solution. 


H. PROPAGATION OVER A PEC CLIFF 
Measured data for the magnitude of the field strength over an aluminum cliff edge was 
also given in [Ref. 10]. Measurements were carried out over the aluminum cliff which was 


situated 10.24/ from the transmitter and of height 1.054 as shown in Fig. 33. 


Transmitter Initial Data Line 
I 





Figure 33. Perfectly conducting cliff edge of height 1.05A. 
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We set up our computation model based on the geometry given in Fig. 33. At a 
frequency = 9.6 GHz, we placed a vertical line source at 10.244 from the edge of the cliff, - 
and set the initial data line at 5.124. We have chosen Ax = Ay = 0.08A (instead of Ax = Ay 
= 0.14) to enable plotting of data at a height, A, = 0.164 (2 vertical steps) at which height the 
measurements were also taken. The result of the numerical solution is plotted in Fig. 34. From 
the plot, we observe that the field attenuates over the first flat section as if the cliff is were 
absent. Immediately beyond the cliff, however, the field strength drops suddenly in the 
shadow region, and then increases with further increase in range, before returning as before 
to decreasing values with increased range. We have computed the attenuation factor based 
on the measured data given in Fig. 6, pp. 39 of [Ref. 10] using Eq. (78) and superimposed 
_ them in Fig. 34. Again, there is excellent agreement between the numerical solution and the 


experimental results, even for points near the cliff edge and immediately beyond the cliff. 


I. PROPAGATION OVER A DIELECTRIC COATED CLIFF 

Experimental results were also available for a dielectric coated cliff [Ref.10]. 
Measurement have been carried out over the same cliff edge given in Fig. 33, but with the 
horizontal surfaces constructed from polypropylene covered aluminum sheet and the vertical 
edge made of aluminum sheet alone. The normalized surface impedance of the polypropylene- 
covered sheet with e, = 2.2 and dielectric thickness, can be calculated to be A = j0.171, 
which is highly inductive. We use the same computational model as in Section H and plot the 
results in Fig. 35. Again, we compute the attenuation factor based on the measured data given 
in Fig. 7, pp. 40 of [Ref. 10] using Eq. (78) and superimposed it in Fig. 35. Again, we observe 
very good agreement between the magnitude of the of the attenuation factor between the 
numerical solution and the measured value. This is true for points near the cliff edge and 
immediately beyond the cliff edge. The behavior of the field strength computed numerically 
also corresponds very well with the experimental results. The results in Sections H and | 
clearly demonstrated the versatility of this numerical method for terrain with PEC surface and 


dielectric coated surface. 
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Figure 34. Magnitude of the normalized fields across a perfectly conducting cliff. 
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Figure 35. Magnitude of the normalized fields across dielectric coated horizontal 
surfaces on a cliff edge; €, = 2.2. Dielectric thickness 0.15 cm. Normalized surface 
impedance is j0.171. 
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J. PROPAGATION OVER A MIXED PATH (LAND-SEA-LAND) 

Next, we would like to examine the capability of the parabolic equation formulation 
to predict radiowave propagation over long range and inhomogeneous ground. Consider a 
four-thirds earth atmosphere over a ground with electrical properties that change along the 
propagation path. The initial and final portions of this "mixed-path" problem are over land (e€, 
= 15, o =0.03 mho/m), while from 47 to 50 km from the source the propagation takes place 
over sea (€, = 80, o = 4 mho/m). Results of path-loss predictions at a receiver height z, = 5 
m for a 59.7 MHz transmitter at a height of z,= 5 m using IFDG method with y, = 10m, x, 


+ ini 


= 40 km, Ay = 1m, and Ax = 100 m for range, r < 47 km and r > 52 km, Ax = 10 m for 47 


km <r < 50 km; and Ax = 1 m for 50 km <r < 52 km is available in [Ref. 1] and the 


Millington approximation in [Ref. 13]. We computed our numerical solution based on the 
geometry described in [Ref. 1]. In view of the large distance, the asymptotic form of Hankel 


function was used in Eq. (76) to compute the magnetic field at the initial range 


HR) = | eX 309, (79) 
| 7™R 


where R = |k.(x-x,)|. The path loss is computed as 


L = -20 log 1% ¢ (4B), (80) 
r 


and plotted in Fig. 36. We see that there is excellent agreement in the path loss (attenuation 
factor) between our computed results and the those computed using the IFDG method [Ref. 
1]. We also noticed that the general behavior of our numerical solution is identical to the 


prediction given in [Ref. 1]. 


K. PROPAGATION OVER A LOSSY VALLEY 
We also investigate the capability of the parabolic equation formulation to predict 
radiowave propagation over a lossy valley. Consider a vertical source radiating on a smooth 


horizontal ground at a distance 30 km from the edge of a 100-m-deep valley that is symmetric 
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about the midpoint of its horizontal bottom, which is 1 km long. The slope of each wall is 
+1/10, so that the total horizontal length of the terrain irregularity is 3 km. The ground is 
characterized by e«, = 4, o = 0.001 mho/m. Results using the IFDG prediction and the 
WAGNER method [Ref. 12] for receiver on the ground are given in [Ref. 1]. In [Ref. 1], the 
IFDG results were obtained using x,,,= 25 km, Ay = 1m, Ax = 20 m for all range, r, except 
for 30 km <r <35 km for which Ax = 1 m. Again, our computational geometry is similar to 
the IFDG geometry. The initial field is computed using Eq. (76) with the asymptotic form for 
Hankel function given in Eq. (79). The path loss for propagation at a frequency of 10 MHz 
is Computed using Eq. (80) and plotted in Fig. 37. Although there is a slight deviation in the 
path loss between our numerical results and the IFDG method [Ref.1], we notice that the 
general trend of the path loss between our method and the IFDG method remains identical. 
The deviation in the magnitude of the path loss could be due to the fact that [Ref. 12] has 
considered an azimuthally symmetrical profile while our profile is 2-dimensionai and laterally 


symmetrical. 
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-~—— Numerical Solution 


O IFDG Method [Ref. 1] 
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Figure 36. Path loss over a land-sea-land path, f= 59.7 MHz, transmitter at 5m, receiver 
at 5m, vertical polarization. 
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Figure 37. Path loss for propagation over a valley, f= 10 MHz, transmitter and receiver 
on the ground. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


In this thesis, a numerically efficient method to model tropospheric radiowave 
propagation over irregular terrain was implemented and tested. The parabolic equation 
method was used for predicting radiowave propagation. A tropospheric boundary condition 
was used to truncate the computational domain at the top, while an impedance boundary 
condition was used on the uneven terrain to characterize the ground. Because the parabolic 
equation method is a full-wave method, it includes all aspects of wave propagation such as 
forward reflection, refraction, diffraction, and surface wave propagation. However, it ignores 
back-scattering. Since the parabolic equation method models wave propagation only in the 
forward direction, it allows for a rapid solution of the fields by way of marching along in 
range starting from an initial range. A great advantage of the parabolic equation method 
compared to the commonly used ray method is that it is valid in the shadow region where the 
latter method completely breaks down. Furthermore, it appears to be the only practical 
method for predicting propagation over long ranges (thousands of wavelengths) and over a 
wideband (HF through SHF). 

The main advancement made in this thesis is the use of a non-rectangular mesh in a 
curvilinear coordinate system to model the uneven terrain. A body fitted coordinate is 
generated based on the specification of the ground profile, and the parabolic PDE, together 
with the boundary conditions, are cast in a curvilinear coordinate system. This not only makes 
the method more efficient, but also permits more accurate imposition and modeling of the 
boundary conditions. We used the Crank-Nicolson implicit scheme in the computational 
domain to solve the parabolic PDE. For a cartesian grid, this method 1s convergent and stable 
for all finite values of step size. In this method, the parabolic PDE 1s considered as being 
satisfied at the mid-point of the computational grid. Central difference approximations, which 
are second-order accurate, are used in the interior as well as at the boundaries. 

Numerical results to predict radiowave propagation over various obstacles were 


computed and validated with the results available in the open literature. These include both 
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a PEC plane and impedance ground with differing permittivity and conductivity. Excellent 
agreement was observed and demonstrated for the following cases : a PEC circular boss, 
lossy Gaussian hill, PEC triangular hill, loaded and unloaded cliff edges, mixed path, and a 
lossy valley. Three-dimensional obstacles are outside the scope of this thesis. 

The model presented in this thesis can be applied to problems in short range 
communications or radar target detection. At HF/VHF frequencies, the antennas are 
electrically close to the ground leading to significant perturbation of the received signals. 
Assessment of propagation of HF/VHF signals over inhomogeneous and irregular terrain is 
highly desirable. In the case of a radar system, signals received from a target depend on the 
direct ray and a ground reflected ray. The path taken by the ground reflected ray depends on 
‘the terrain roughness and ground constants. Since target detection is based on the composite 
signal, it is important to assess ground conditions and atmospheric factors. The numerical 
model developed in this thesis provides a useful tool for predicting the performance of radar 
or communication links before commencing to design, develop, and deploy the system in the 
field. From the model developed, one can predict propagation factors such as path loss, 
coverage diagrams, and the perturbation of antenna radiation patterns over realistic terrains 
and ranges. 

This thesis represents the beginning of an effort to predict radiowave propagation 
using the parabolic equation method. In its standard form, the accuracy of the method is 
limited to waves traveling within +15° from the horizontal. It is recommended that follow-on 
work should use a wide-angle parabolic equation method to accommodate waves traveling 


over wider angles and in an inhomogeneous atmosphere. 
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APPENDIX A. MATLAB SOURCE LISTINGS 


The formulation described in Chapter II and the numerical procedure described in 
Chapter IJ] were implemented using MATLAB sofware. The routines that generated the 


results in Chapter IV are listed below. 


A. EXECUTIVE ROUTINE 


% Filename : main.m 

% Title : Main Program or Executive Program 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% This is the main program for computing the Parabolic 
% Equation Solution of Radiowave Propagation in an 

% Inhomogenous Atmosphere and Over Irregular Terrain. 
% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% U : computed Hz-Field 

% wavelength : wavelength (frequency/3e8) 

% time : total CPU time used 

% 

% Associated Matlab Files : 

% Called by none 

% Subroutines : data 

% dinp 

“% cmet 

% hfld ‘ 

% epde 

% 

% Associated Functions =: none 

% 

clear all, 

t=cputime; % start CPU clock 

data % input data file 

dinp % compute the physical domain 
cmet % compute the metric 

hfld % compute Hz-Field on initial data line 
epde % march and compute Hz-Field 


time=cputime-t: 
save output xphy yphy N M U time wavelength, 
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B. PHYSICAL DOMAIN GRID 


% Filename : dinp.m 

“ Title : Computation of Physical Domain 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 


% Comments : 
% This subroutine computes the entire computational 
% domain based on the user input data for the physical 
% ground plane. | 
% 
“ Input Variables : 
% xin : x-coordinate of ground plane 
% yin : y-coordinate of ground plane 
% 
% Output Variables : 
% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% 
% Associated Matlab Files : 
% Called by > main 
% Subroutines : none 
% 
% Associated Functions _: none 
% 
nseg=length(xin)-1; % total number of segment 
% 
%%%7%o%. Compute horizontal axis of the physical domain %%»%%% 
% 
st=1; 
for 1=] :nseg: 
npts(1)=round(sqrt((xin(i+1 )-xin(4))2+(yinGt1 )-yin(1))*2)/xdel); 
xphy(st)=xin(i); 
yphy(st)=yin(i); 


for j=] :npts(1)-1; 
xphy(st+j)=xin(i)+j*(xin(i+1)-xin(i))/npts(i); 
yphy(1,st+j)=yin(i)+j*(yin(i+1)-yin(i))/mpts(i): 
end,  - 
st=st+npts(1); 
end: 
xphy(st)=xin(length(xin)); % include the last point for x 
yphy(st)=yin(length(yin)), % include the last point for y 
% 
AAAY%e Compute vertical axis of the physical domain %%%%% 
% 
M=length(xphy)- 1; 
N=round(yup/ydel); 
for 1=1:N; 
yphy(+1,:)=(1 -/N).*yphy(1,:)+@/N)*yup; 
end; 
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C. CURVILINEAR COORDINATE SYSTEM METRIC 


% Filename : cmet.m 

% Title : Generation of Curvilinear Coordinate System Metric 
% Revision No.:RI.0 _ 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% This subroutine computes the various metrics or derivatives 
% (dx_x1, dx_eta, dy_xi, dy_eta, ddy_eta) at any interior 
% points. In the Crank Nicolson implicit scheme, the metrics 
% are needed at the midpoint with respect to xi, 1.¢., at 
% xi=xitdel_xi/2. 

% 

% Input Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% 

% Output Variables : 

% dx_xi : Ist order x-derivative w.r.t. xi 

% dx_eta : Ist order x-derivative w.r.t. eta 

% dy_xi : Ist order y-derivative w.r.t. x1 

% dy_eta : Ist order y-derivative w.r.t. eta 

% ddy_eta : 2nd order y-derivative w.r.t. eta 
% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 

% 

% Associated Functions =: none 

% 


msize=ones(size(yphy(:,1))); 
dx_xi=msize*(xphy(2:M+] )-xphy(1:M)); 
dx_eta=0; 
for 1=1:N+1,; 
dy_xi(i,]:M)=(1-G-1)/N).*(yphyd,2:M+1 )-yphy(, 1:M)); 
end; 
dy_eta=msize*((yphy(N+1,1)-(yphy(1,2:M+1)+yphy(1,1:M))./2)./N); 
ddy_eta=0; 


D. MAGNETIC FIELD ON INITIAL DATA LINE 


% Filename : hfld.m 

% Title : Generation of Hz-Field on Initial Data Line 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% This subroutine generates the magneic field (Hz-Field) 
% at the initial data line for a line source over 

% an impedance plane. 

% (Note : Never place the initial data line at the 
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% same location as the transmitter, else the field 


% generated would be zero). 

% 

% Input Variables : 

% yphy : y-coordinate of physical domain 
% 

“% Output Variables : 

% Hz © : Magnetic Field @ initial data line 
% 

% Associated Matlab Files : 

% Called by > main 

% Subroutines : none 

% 

% Associated Functions : hankel2 

% 


yinit=[yphy(1:N,1)' yphy(N+1,1):ydel:ytop]; 

Ro=sqrt((xphy(1 )-x0).%2+(yinit-y0).%2)*ko; 

% 

WAVY Generate Hz-Field using Fredholm Equation %%%%% 
% 


hr2=hankel2(0,Ro); % load Hankel function of 2nd order 
H=2*((1j}*ko)/4).*((xphy(1)-x0)./Ro).*hr2; 

Hz=H."; 

% 

“ Note: "." is for Nonconjugated Transpose. Using only "" 

% will result in conjugated transpose (i.e. the sign 

% of the complex term is inverted) 


E. FORMULATION OF PARABOLIC EQUATION MATRIX 


% Filename : epde.m 

% Title : Matrix Equation for the Parabolic System Equations 
% Revision No. : R1.0 

‘% Date of Last Revision : 10 Mar 95 

% Comments : 


% This program defines the PDE for the matrix equations 
% and performs the numerical implementation using the 
% Crank-Nicolson Implicit Scheme (Discretization). 

% The applicability of the PDE equation is also 

% extended to accommodate the derivative boundary 

% conditions on the lower and upper boundaries. 

% 

% Input Variables : 

% dx_xi : lst order x-derivative w.r.t. xi 

% dx_eta : Ist order x-derivative w.r.t. eta 

% dy_xi! : Ist order y-derivative w.r.t. xi 

% dy_eta : Ist order y-derivative w.r.t. eta 

% ddy_ eta : 2nd order y-derivative w.r.t. eta 
% Hz : Magnetic Field @ initial data line 
% 
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% Output Variables : 


% U : computed Hz-Field 
% 

% Associated Matlab Files : 

% Called by > main 

% Subroutines - elbc 

% eubc 

% 

% Associated Functions _: none 

% 


AW%V%%% Assume constant atmospheric refractive index %*%%%% 

n=]; dn_x=0; dn_y=0; 

al=-(2/n)*dn_x; a2=-(2/n)*dn_y; 

al_star=(n‘*2-]-2)*(al/ko))*(ko%2),; 

% 

%%%%% Set up the matrix coefficients, (M x N) for b1, b2, b3 %*%%%% 
% 

temp 1=2)*ko-al ; 

-b1]=dx_xi.*al_star./temp 1, 
b2=dx_x1./dy_eta.*((a2-ddy_eta./dy_eta.*2)./temp1+dy_x1./dx_x1); 
b3=dx_x1./(temp!.*dy_eta.%2); 
alpha=0.5.*(b3+0.5.*b2); 
beta=(b3-0.5.*b1); 
gamma=0.5.*(b3-0.5.*b2); 

% 

%%%%% Initialization of Magnetic Field %%%%% 

% 

UC1:N+1,1I=Hz(1:N+1); 

% 

%%%%% Compute coefficient of Lower BC %%%%% 
% 

elbc 

% 

%%%%% Compute invariant coefficient of Upper BC %%%%% 
% 

r(1:M)=sqrt((8j*ko)/pi)./sqrt((xphy(2:M+1 )-xphy(1:M))./2); 
% 


%%%%% Loading the matrix %%%%% 

% 
for p=1:M; 

A=spdiags([gamma(:,p) -(1+beta(:,p)) alpha(:,p)], -1:1, N+1, N+1); 
B=spdiags([-gamma(:,p) -(1-beta(:,p)) -alpha(:,p)], -1:1, N+1, N+1); » 
A(1,2)=alpha_p(1,p); A(I,1)=-(1+beta_p(1,p)); 
B(1,2)=-alpha_p(1,p); B(1,1)=-(1 -beta_pp(1,p)): 

% 

%%%%% Compute coefficient of Upper BC %**%*%%% 

% 

eubc, 

lambda=beta(N+1,p)+dy_eta(p)*2*alpha(N+1,p)*(r(p)); 

gamma _ p=gamma(N+1,p)t+alpha(N+1,p); 
A(N+1,N+1)=!+lambda; 

A(N+1,N)=-gamma_p; 
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BC(N+1,N+1)=(1 -lambda); 
B(N+1,N)=gamma_p; 

B(N+1,N+2)=1; 
U(N+2,p)=4*alpha(N+1,p)*dy_eta(p)*(s(p)); 
U(1:N+1,p+1)=(A\B)*U(1:N+2,p); 

end; 


F. LOWER BOUNDARY CONDITION 


% Filename : elbc.m 

% Title : Lower Boundary Conditions 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 
% Comments : 


% This program formulates the equations for the 

% impedance boundary conditions imposed on the irregular 
% lower boundary. Impedance boundary conditions relates 
% the tangential components of electric and magnetic 

% fields at the interface of two media. 

% The surface impedance is determined from the 

% intrinsic impedance of the medium by considering 

% plane wave reflections from the interface. 

% 

% Input Variables : 

% dx_xi : Ist order x-derivative w.r.t. xi 

% dy_ xi : Ist order y-derivative w.r.t. xi 

% dy_eta : Ist order y-derivative w.r.t. eta 

% alpha : coefficient for U(p,qt1) term 

% gamma : coefficient for U(p,q-1) term 

% 

% Output Variables : 

% alpha_p : coefficient for U(p,1) term 

% beta_p : coefficient "(1+beta_P)" for U(p,0) term 
% beta_pp : coefficient "(1-beta_pp)" for U(p-1,0) term 
% 

% Associated Matlab Files : 

% Called by : epde 

% Subroutines : none 

% 

% Associated Functions —_: none 

% 


temp2=sqrt(dx_xi(1,1:M).42+dy_xi(1,1:M).%2); 
ct=dx_xi(1,1:M)./Atemp2; % cos(theta) 
st=dy_x1i(1,1:M)./temp2; % sin(theta) 
c]l=-1)*ko*(n“2.*delta_v - st); 
alpha_p=alpha(1,:)+gamma(1,:); 

temp3=2.* gamma(1,:).*dy_eta(1,:).*ct; 
beta_p=beta(1,:)-temp3.*(cl -2*st./dx_xi(1,1:M)); 
beta_pp=beta(1,:)-temp3.*(c1+2*st./dx_xi(1,1:M)); 
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G. UPPER BOUNDARY CONDITION 


% Filename : eubc.m 

% Title : Upper Boundary Conditions 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 
% Comments : 


% This program formulates the equations for the upper 

% boundary conditions at the tropospheric boundary. 

% To truncate the computational domain at the upper 

% boundary, a point high enough where the atmosphere 

% is homogeneous is considered. 

% : 

% Input Variables : 

% dx_xi : 1st order x-derivative w.r.t. x1 

% dy_ x1 : Ist order y-derivative w.r.t. xi 

% dy_eta : 1st order y-derivative w.r.t. eta 
% alpha : coefficient for U(p,qt+1) term 

% gamma : coefficient for U(p,q-!) term 

% 

% Output Variables : 

% s : coefficient for U(p,1) term 

% 

% Associated Matlab Files : 

% Called by : epde 

% Subroutines : none 

% 

% Associated Functions _ : fresnel 

% 


tempA=0; tempB=0; tempC=0; 

x=(xphy(p+1)+xphy(p))/2; % x(p-1/2) 

xd=(xphy(p+1)-xphy(p))/2; % x(p-1/2)-x(p-1) 

for 1=N+2:round(ytop/ydel+1); 
j=i-N-1; 
temp4(j)=fresnel(sqrt(ko/(pi*(x-xin(1)))).*(ymit(i)-yimit(N+1))); 
temp5(j)=fresnel(sqrt(ko/(p1*(x-xin(1)))).*(yimit@- 1)-yimit(N+1))); 
tempA=tempA+(Hz(i,1)-Hz(i-1,1))/ydel.*(temp4(j)-temp5(j)); 

end; 

% 

%%%%% Compute the UCN+1,p-1) term in s(p) %%%%% 

% 

tempA=sqrt(2})*tempA, 

tempC=sqrt(8j*ko/pi)*U(N+1 ,p)/sqrt(xd), 

if p=1, 
s(p)=tempAttempC; 

else; 
temp6=sqrt(x-xphy(2:p))+sqrt(x-xphy(1 :p-1)), 
tempB=sum((U(N+1,2:p)-U(N+1, 1 :p-1))./temp6), 
s(p)=tempA-sqrt(8j*ko/pi)*tempB+tempC, 

end, 
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H. FRESNEL INTEGRAL FUNCTION 


% Filename : fresnel.m 

% Title : The Fresnel Integral Function 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 
% Comments : 

% 

% Input Variables : 

% 

% Output Variables : 

% 

% Associated Matlab Files : 

% 

% Associated Functions : 

% 

function F=fresnel(limits); 


F=0,; 
elseif limits < 0; 

limits=- 1 *limits; 
tau=0:interval:limits; 

tau=-1].*tau; 
integral=exp(-1j*pi/2*(tau.*tau)); 
F=trapz(tau,integral); 
else 

tau=O:interval:limits; 
integral=exp(-1)j*pi/2*(tau.*tau)); 
F=trapz(tau,integral); 
end; 


I. INPUT DATA FILE FOR PROPAGATION OVER PEC PLANE 


% Filename : pec.m 

% Title : Input Data File 

% Revision No. : R1I.0 

% Date of Last Revision : 10 Mar 95 
% Comments : 


% This is the input data file for computing the Parabolic 
% Equation Solution of Radiowave Propagation in an 

% Inhomogenous Atmosphere and Over Irregular Terrain. 
% It contains the frequency of operation, the x and y 

% discretization size, x- and y-coordinates of the ground 
% plane domain, and the upper boundary. 

% It also contains the surface impedance (Note : surface 
% impedance = 0 for PEC). 

% 

% Input Variables : none 

% 
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% Output Variables : 


% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 
% 
% Associated Matlab Files : 
% Called by : main 
% Subroutines : none 
% 
% Associated Functions =: none 
% 
freq=300e6, % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0. 1 * wavelength; % size of delta-x 
ydel=0. 1 *wavelength; % size of delta-y 
% 
%%wH%%% Define Ground Plane Coordinates %*%%%% 

Y 


xin(1)=5*wavelength; yin(1)=0; 
xin(2)=15*wavelength; yin(2)=0; 


yup=20*wavelength; % y-upper boundary 
ytop=5*wavelengthtyup; % troposphere boundary (@ infinity) 
x0=0; yO=0; % position of transmitter 

% 

%U%%%% Define Surface Impedance %%%%% 

% 

delta_v=0; % surface impedance (PEC) 


J. INPUT DATA FILE FOR PROPAGATION OVER PEC KNIFE-EDGE 


% Filename : kpec.m 

% Title : Input Data File for PEC Knife Edge 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 

% 

% Associated Functions none 

% 
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freq=300e6; % frequency of operation 


wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0. 1 *wavelength; % size of delta-x 
ydel=0.1*wavelength; % size of delta-y 
% 

AVY Define Ground Plane Coordinates %*%%%% 
% 


xin(1)=5*wavelength; yin(1)=0, 
xin(2)=9.9*wavelength; yin(2)=0; 
xin(3)=10*wavelength; yin(3)=2*wavelength: 
xin(4)=10.1*wavelength; yin(4)=2*wavelength: 
xin(5)=10.2*wavelength;yin(5)=0; 
xin(6)=30*wavelength; yin(6)=0; 


yup=22* wavelength; Yo y-upper boundary 
ytop=5*wavelengtht+yup; % troposphere boundary (@ infinity) 
x0=0; yO=0; % position of transmitter 

% 

WAVY Define Surface Impedance %%%%% 

% 

delta_v=0; ‘% surface impedance (PEC) 


K. INPUT DATA FILE FOR PROPAGATION OVER CIRCULAR BOSS IN A 
PEC PLANE 


% Filename : boss.m 

% Title : Input Data File for Semi-Circular Boss over PEC Plane 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : maim 

% Subroutines : none 

% 

% Associated Functions _: none 

% 

freq=3e6, % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0.0 1 *wavelength; % size of delta-x 
ydel=0.05*wavelength; % size of delta-y 

% 
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%%%%% Define Ground Plane Coordinates %%%%% 

% 

b=0.5*wavelength; 

A=wavelength, B=0.8*wavelength; 
x_mid=-0.5*wavelength:0.01*wavelength:0.5*wavelength; 
y_mid=sqrt(b*2-x_mid.“2), 

x_fnt=-0.6*wavelength:0.01 *wavelength:-0.5 1 *wavelength; 
y_fnt=zeros(size(x_fnt)); 

x_bk=0.51*wavelength:0.01 *wavelength:wavelength; 
y_bk=zeros(size(x_bk)); 

xin=[x_fnt x_mid x_bk]; 

%xin=xin+0.75*wavelength; 

yin=[y_fnt y_mid y_bk]; 

clear x_fnt; clear x_mid; clear x_bk; 

clear y_fnt; clear y_mid; clear y_bk; 


yup=20* wavelength; % y-upper boundary 
ytop=5*wavelength+yup; % troposphere boundary (@) infinity) 
x0=-0.75*wavelength; yO=0.1*wavelength; % position of transmitter 

% 

%%%%% Define Surface Impedance %*%%/% 

% 

delta_v=0; % surface impedance (PEC) 


L. INPUT DATA FILE FOR PROPAGATION OVER LOSSY IMPEDANCE 
PLANE 


% Filename : imped.m 

% Title : Input Data File for Lossy Impedance Plane 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 

% 

% Associated Functions =: none 

% 

freq=1e6, % frequency of operation 
wavelength=3e8/freq; _ % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0. 1 *wavelength; % size of delta-x 
ydel=0. | *wavelength, % size of delta-y 
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% 

%%%%% Define Ground Plane Coordinates %*%%%% 
% 

xin(1)=0.5*wavelength; yin(1)=0; 

xin(2)=20* wavelength; yin(2)=0: 


yup=20* wavelength; % y-upper boundary 
ytop=5*wavelengthtyup, % troposphere boundary (@ infinity) 
x0=0; yO=wavelength/300,; % position of transmitter 

% 

WAVYY Define Surface Impedance %%%%% 

% 

er=10; % relative permittivity 
sigma=10e-3;, % conductivity 


sigma_1=60*wavelength*sigma, § % relative conductivity 
delta_v=sqrt(1/(er-1j*sigma_r));% surface impedance (vertical pol) 


_M. — INPUT DATA FILE FOR PROPAGATION OVER LOSSY GAUSSIAN HILL 


“% Filename : gauss.m 

% Title : Input Data File for Lossy Gaussian Hill 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% 

“% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 

% 

% Associated Functions =: none 

% 

freq=1e6, % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0. 1/3* wavelength; % size of delta-x 
ydel=0. 1 *wavelength; % size of delta-y 

% 

%M%%%Y% Define Ground Plane Coordinates %%%%% 

% 

Ht=1000; c=5000; w=3000, % constants for gaussian hill 


xin=2000:xdel:24*wavelength;% x-coordinate 

yin=Ht. *exp(-9.*(((xin-c)./Aw).*2));, % y-coordinate 

yup=30* wavelength; % y-upper boundary 
ytop=5*wavelengthtyup; % troposphere boundary ((@) infinity) 
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x0=0; yO=wavelength/100; % position of transmitter 


% | 

%%%%% Define Surface Impedance %%%%% 

% 

er=10; % relative permittivity 
sigma=10e-3; % conductivity 


sigma _r=60*wavelength*sigma, §% relative conductivity 
delta_v=sqrt(1/(er-1)*sigma_r));% surface impedance (vertical pol) 


N. INPUT DATA FILE FOR PROPAGATION OVER PEC ISOSCELES HILL 


% Filename : 1sos.m 

% Title : Input Data File for PEC Isosceles Hill 

% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

“% Comments : 

% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 
% 

% Associated Matlab Files : 

% Called by > main 

% Subroutines : none 

% 

% Associated Functions _: none 

% 

freq=9.6e9; % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0.08* wavelength; % size of delta-x 
ydel=0.08* wavelength; % size of delta-y 

% 

%%%%% Define Ground Plane Coordinates %*%%%% 
% 

xin(1)=5.12*wavelength, yin(1)=0; 
xin(2)=10.24*wavelength; yin(2)=0; 

xin(3 )=20.32* wavelength; yin(3)=1.232* wavelength; 
xin(4)=30.4* wavelength, yin(4)=0; 
xin(5)=64*wavelength; yin(5)=0; 


yup=20* wavelength; % y-upper boundary 
ytop=5*wavelengthtyup; % troposphere boundary (@ infinity) 
x0=0; yO=0; % position of transmitter 

% 

%w%V%%%% Define Surface Impedance %%%%% 

% 

delta_v=0; % surface impedance (PEC) 
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O. INPUT DATA FILE FOR PROPAGATION OVER PEC CLIFF 


% Filename : cliff.m 

% Title : Input Data File for PEC Cliff 

% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 

% 

‘% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 
% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 

% 

% Associated Functions _: none 

% 

freq=9.6e9; % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0.08* wavelength; % size of delta-x 
ydel=0.08* wavelength; % size of delta-y 

% 

AVY%Y%% Define Ground Plane Coordinates %%%%% 
% 

xin(1)=5.12*wavelength; yin(])=1.0464*wavelength; 
xin(2)=10.24*wavelength; yin(2)=1.0464*wavelength; 
xin(3)=10.34*wavelength; yin(3)=0: 

xin(4)=80* wavelength; yin(4)=0; 


yup=2 | *wavelength, Y% y-upper boundary 
ytop=5*wavelength+yup; ‘% troposphere boundary (@ infinity) 
x0=0; yO=yin(1); % position of transmitter 

% 

WVY%V%Y% Define Surface Impedance %%%%% 

% 

delta_v=0; % surface impedance (PEC) 


P. INPUT DATA FILE FOR PROPAGATION OVER DIELECTRIC COATED 
CLIFF 


% Filename : diele.m 

% Title : Input Data File for Dielectric coated Cliff 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 
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% 
% Input Variables : none 


% 

% Output Vanables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : None 

% 

% Associated Functions : none 

% 

freq=9.6e9, % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
xdel=0.08* wavelength; % size of delta-x 
-ydel=0.08* wavelength; % size of delta-y 


% 

%%%%% Define Ground Plane Coordinates %%%%% 
% 

xin(1)=0.5*wavelength, yin(])=1.0464*wavelength; 
xin(2)=10.24*wavelength; yin(2)=1.0464*wavelength; 
xin(3)=10.34*wavelength; yin(3)=0; 

xin(4)=80* wavelength; yin(4)=0; 


yup=2 | *wavelength; % y-upper boundary 
ytop=5*wavelength+yup; % troposphere boundary (@ infinity) 
x0=0; yO=yin(1); % position of transmitter 

% 


%%%%% Define Surface Impedance %*%%%% 
% 
delta_v=0.171); % surface mpedance 


Q. INPUT DATA FILE FOR PROPAGATION OVER MIXED PATH 


% Filename : mix.m 

% Title : Input Data File for Land-Sea-Land 
% Revision No. : R1.0 

% Date of Last Revision : 10 Mar 95 

% Comments : 


% 

% Input Variables : none 

% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 


% Associated Matlab Files : 
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% Called by : main 

% Subroutines : none 

% 

% Associated Functions _: none 

% 

freq=59.7e6; 

wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; Yo Free space wave number 

ydel=0.2; % y-delta size 

% 

AAV%% Define Ground Plane Coordinates %%%%% 

% 

xin(1)=40000; yin(1)=0; xdel(1)=100, impz(1)=sqrt(1/(15-9.045j)); 
xin(2)=47000; yin(2)=0; xdel(2)=10: impz(2)=sqrt(1/(80-1206))); 
xin(3)=50000; yin(3)=0, xdel(3)=1: , impz(3)=sqrt(1/(15-9.045j)); 
xin(4)=52000, yin(4)=0, xdel(4)=100; impz(4)=sqrt(1/(15-9. 045))), 
xin(5)=60000; yin(S)=0; impz(5)=sqrt(1/(15-9.045j)): 

yup=10; 

_ytop=5*wavelengthtyup; 

x0=0; yO=5; 


R. INPUT DATA FILE FOR PROPAGATION OVER LOSSY VALLEY 


“% Filename : valley.m 

% Title : Input Data File for Valley 
% Revision No. : RO.2 

% Date of Last Revision : 20 Feb 95 
% Comments : 

% 

Input Variables : none 


% 

% Output Variables : 

% xphy : x-coordinate of physical domain 
% yphy : y-coordinate of physical domain 
% wavelength : wavelength (frequency/3e8) 

% 

% Associated Matlab Files : 

% Called by : main 

% Subroutines : none 


% 
% Associated Functions : none 


% 

freq=10e6; % frequency of operation 
wavelength=3e8/freq; % wavelength 
ko=2*pi*freq/3e8; % free space wave number 
ydel=0. 1 *wavelength, % size of delta-y 

% 

Y%V%%% Define Ground Plane Coordinates %*%%%% 

% 


xin(1)=28000; yin(1)=100; xdel(1)=300; 
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xin(2)=30000; yin(2)=100; xdel(2)=3; 
xin(3)=3 1000; yin(3)=0; xdel(3)=3; 
xin(4)=32000; yin(4)=0; xdel(4)=3; 
xin(5)=33000; yin(5)=100, xdel(5)=3;, 
xin(6)=35000; yin(6)=100; xdel(6)=30; 
xin(7)=50000; yin(7)=100; 
yup=50+yin(1); 
ytop=5*wavelengthtyup; 


x0=0, yO=yin(1); 

% 

WAVA%Y% Define Surface Impedance %%%%% 

% 

er=4, % relative permittivity 
sigma=1e-3; % conductivity 


sigma r=60*wavelength*sigma, §% relative conductivity 
delta_v=sqrt(1/(er-1j*sigma_r));% surface impedance (vertical pol) 
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